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Abstract 


A meshless local Petrov-Galerkin (MLPG) method that uses radial 
basis functions (RBFs) as trial functions in the study of Euler-Bernoulli 
beam problems is presented. RBFs, rather than generalized moving least 
squares (GMLS) interpolations, are used to develop the trial functions. 

This choice yields a computationally simpler method as fewer matrix 
inversions and multiplications are required. Test functions are chosen 
as simple weight functions as they are in the conventional MLPG 
method. Compactly and noncompactly supported RBFs are considered. 

Noncompactly supported cubic RBFs are found to be preferable. Patch 
tests, mixed boundary value problems, and problems with complex 
loading conditions are considered. Results obtained from the radial 
basis MLPG method are either of comparable or better accuracy than 
those obtained when using the conventional MLPG method. 

Introduction 

Meshless methods are developed to overcome some of the disadvantages of the finite element method 
(FEM), such as discontinuous secondary variables across interelement boundaries and the need for 
remeshing in large deformation problems (Atluri, Cho, and Kim, 1999; Atluri and Shen, 2002b; Atluri 
and Zhu, 1998; Belytschko, Liu, and Gu, 1994; and Nayroles, Touzot, and Villon, 1992). Recent litera- 
ture shows extensive research on meshless methods and, in particular, the meshless local Petrov-Galerkin 
(MLPG) method. Two recent monographs (Atluri and Shen, 2002a; and Liu, 2002) summarize the status 
to date. Most literature published to date on the MLPG method presents variations of the method for C° 
problems. A comparatively limited amount of work (Atluri, Cho, and Kim, 1999; Donning and Liu, 
1998; Gu and Liu, 2001a and 2001b; and Krysl and Belytschko, 1995) is reported on C 1 problems. The 
additional degrees of freedom involved in C 1 problems introduce levels of complexity that are difficult to 
implement in meshless methods. Gu and Liu present a local point interpolation method (LPIM) for static 
and dynamic analysis of thin beams (Gu and Liu, 2001a). Atluri, Cho, and Kim (1999) present an analy- 
sis of thin-beam problems by using a Galerkin implementation of the MLPG method, in which a general- 
ized moving least squares (GMLS) approximation is used to construct the trial functions, and the test 
functions are chosen from the same space. In Phillips and Raju (2002) and Raju and Phillips (2002, 
2003a, and 2003c), an MLPG method is presented in which a GMLS approximation is used to construct 
the trial functions; test functions are chosen from a different space. Closer scrutiny of these formulations 
shows that a large number of calculations are required to compute the first- and second-order derivatives 
of the moving least squares (MLS) trial functions. Flence, a computationally simpler alternative to the 
MLS trial functions that provides the same accuracy as the MLS functions is preferred. 

In this paper, the use of radial basis interpolation functions (see Powell, 1992; Wendland, 1999; and 
Wu, 1995) as trial functions is explored in the MLPG formulation for beam problems. The radial basis 
functions (RBFs) and their implementation are simple, and the evaluation of the derivatives is much 
simpler than for the traditional MLS approximations (Raju, Phillips, and Krishnamurthy, 2003b). In the 
present radial basis MLPG formulation (abbreviated hereafter as RPG), simple weight functions are 
chosen as test functions, and Gaussian quadrature is used to integrate the weak form. The effectiveness of 
the RPG method is evaluated by applying the formulation to a variety of patch test and mixed boundary 
value problems. 



The intent of this paper is to present an alternative to the MLS interpolation that may be more compu- 
tationally efficient in the MLPG method for beam problems. Currently, the MLPG method is more 
computationally expensive compared to the FEM. The objective of this work is not to improve the 
computational efficiency of the MLPG method but rather to explore the feasibility of using RBFs as 
interpolations in the method. 

The outline of the paper is as follows: First, the MLS interpolation used in the conventional MLPG 
method is discussed as motivation for finding a simpler alternative. Next, an overview of RBFs for C° 
problems is presented, the shape functions obtained from radial basis interpolation are derived, and the 
shape functions that are obtained when polynomial basis functions are included in the interpolation are 
derived. The MLS and radial basis shape function discussions are then expanded and repeated for beam 
(C 1 ) problems. The system of algebraic equations developed from the local weak form of the governing 
differential equation and the chosen trial and test functions is presented. Patch test problems are used to 
validate the RPG method for different choices of RBF, and the RPG method is applied to mixed boundary 
value problems. Finally, the method is applied to problems with complex loading and boundary 
conditions. 

C° Interpolation Schemes 

In this section, the MLS interpolation scheme used in the conventional MLPG method is discussed 
first; then, two interpolation schemes involving RBFs are presented. In the first scheme, RBFs alone are 
used to construct the shape functions. The second scheme is a hybrid that uses both RBFs and polyno- 
mial basis functions to construct the shape functions. 

The MLS Interpolation 

An MLS interpolation is a scheme that fits a smooth function through an assumed set of fictitious 
nodal values. The interpolation is performed such that the sum of the least squares error between the 
interpolated function and the fictitious nodal values is a minimum (Belytschko, Lu, and Gu, 1994; and 
Nayroles, Touzot, and Villon, 1992). A schematic of the MLS interpolation is presented in figure 1 . 



Figure 1. Moving least squares (MLS) interpolation. 


2 


In one-dimensional problems, the C° MLS shape function for node j can be constructed as follows. In 
the equations that follow, n is the number of nodes in the domain of definition of node j, and m is the 
order of the polynomial basis function (Nayroles, Touzot, and Villon, 1992). The shape function (V at 

node j can be written as 



(1) 


where 


(x) is a polynomial basis function 


[p T (x)] = [pi(x), p 2 (x), p 3 (x), ..., p m (x)\ 


I i 2 m — 1 

=1, X, X , ..., X 


( 2 ) 


and 


[A] = [P] T [^][P] 

[B] = [P] T M 


(3) 


In equations (3), [P] is an (n, m ) matrix, and [A,] is a diagonal (n, n) matrix defined as 

-.T 


[P] = 


P T Ul) P T (^2) ••• P T U«) 


Xl(x) 


X = \ 


X 2 {x) 


X n (x) 


(4) 


(5) 


where 'kj(x), j = 1, ..., n, is a weight function. The first derivative of these shape functions is required 
by the MLPG method and is given by (Atluri and Zhu, 1998; and Belytschko, Liu, and Gu, 1994) 


<bj,x 



where ( ) s d ( ) Idx , and 


( 6 ) 


[A ]” 1 =-[A] _ 1 [A] >JC [A] 1 


(7) 


RBF Interpolation 

The radial basis formulation for C° problems is presented first. The RBF provides a continuous inter- 
polating function for «(x) as a linear combination of radial functions (Powell, 1992). The interpolating 
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function is 


TV 

u(x)= ^ Rj(x) cij (8) 

j = 1 

where Rj(x), the RBFs, are functions that are “centered” at each of the TV scattered points, and aj are the 
unknown coefficients, j = 1, 2, . . TV. The RBF, Rj(x), are functions of distance rj and are defined as 

Rj(x)=R j (r J ) (9) 


In two dimensions, the radial distance rj can be expressed in terms of the Cartesian coordinates as 



( 10 ) 


where xj and yj are the coordinates of node j. Equation (8) can be written in matrix form for C° problems 
as 


u(x) = R T (x) a 


(ID 


where 


R T (x) = [R 1 (x), R 2 (x), R 3 (x), ..., R n (x )] 
a = {ci \ , a 2 , a 3 , ..., fljv} T 


( 12 ) 


Forcing the interpolation of equation (8) to pass through the N scattered points, a set of equations to 
determine the coefficients aj in terms of nodal values uj can be written as 


where 




R B a = u 


(13) 

u 

T = («t. 

^2’ 

•» U N ) 

(14) 


" R\ Oh) 

R 2 ( x i) '• 

• %(x l)“ 


R b = 

R\ (x 2 ) 

R 2 (x 2 ) ■■ 

• R n (x 2 ) 

(15) 



R 2 ( x n) ■■ 

■ r n( x n)_ 



Note that R u is an (TV, TV) matrix. Here, u (uj,j = 1, 2, ..., TV) are the nodal values of u at the TV scattered 
points. The unknown coefficients in equation (13) can be obtained as 

a = Rg u (16) 
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The interpolating function for u(\) in equation (11) can now be rewritten as 


N 

M(x) = R T (x)Rg 1 u= ^<pj(x)uj 
j = 1 


The nodal shape functions are then 

cp(x) = R T (x)Ri 1 

= [(pt(x), (p 2 ( x X (p3(x)> 9 aKx)] 


( 17 ) 


(18) 


or 


N 

cp /(x)= X^(x)^/ 


k=l 


(19) 


where are the elements of the matrix . The shape function (py(x) obtained through the above pro- 
cedure satisfies the Kronecker Delta property at the nodes (Atluri and Shen, 2002a), that is, 

(Pj(x k ) = 8jk (20) 

Note that the shape functions in equation (19) also satisfy the property 

N 

Xcp,.(x) = l (21) 

7=1 

only at the nodes. As the number of nodes is increased, however, the shape functions in equation (19) 
satisfy (Raju, Phillips, and Krishnamurthy, 2003b) 

N 

lim X ( P/(x) = 1 (22) 

N^>°° 


everywhere in the domain. 

In contrast to the MLS method, the derivatives of the shape functions are easy to evaluate by using 
equation (19) as 


>/(x) 

dx 


y Mk 

tx * 


t ki 


5ip/(x) 

dy 


V dR k 

It 


tkl 


(23) 
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where 


with 


()R/ C _ dR k dr k 
()x di'j. c)x 

dr k _x-x k 
dx r k 


()Rj c _ dR k dr k 

dy r h dy 


<h _ y~yk 

dy r k 


(24) 


(25) 


Some of the well-known classical radial functions are presented in table 1, in which (tj = fj/sj) and sj is 
some normalizing distance. Note that the shape parameter c in the RBFs (table 1) is user defined. The 
classical radial functions have two limitations: (1) the matrix Rb in equation (15) may not be positive 
definite, and (2) the functions do not possess local support, that is, the changes in a single location (xj, yj) 
affect the entire interpolation. To overcome these limitations, compactly supported positive definite 
radial functions were proposed (Wu, 1995). These functions were derived by using a constraint to 
guarantee positive definiteness of the interpolation matrix Rb. 

Table 1. Classical Radial Basis Functions (RBFs) 


Classical RBF 

Equation 

Linear 

Rj(x) = ct j 

Cubic 

R J (x) = (tj+c ) 3 

Thin plate spline 

R j( x ) = t j log (ctj) 

Gaussian 

Rj(x) = e ct J 

Multiquadric 

Rj (x) = -\Jtj + c 2 


*See also tables 5.2 and 5.3 (Liu, 2002). 


The compactly supported functions assign zero weight to all points not belonging to the compactly 
supported region, thus leading to a sparse interpolation matrix. Adapted compact RBFs are used here: 

Compact-I: 


(1 - tj ) 5 (8 + 40r j + 48 tj + 25 tj + 5 tj), 0 < tj < 1 
0 


(26) 


Compact-II: 


R jW = \ 


(1 - tj ) 6 (6 + 36 tj + 82 tj + lit) + 30F- + 5 tj ), 0 < < 1 

0 t .i >] 


(27) 
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where (tj = rj/sj) and sj, the normalizing distance, is the radius of the domain of compact support. The 
shape functions (eq. (19)) obtained from the Compact-I and Compact-II functions possess all the 
properties in equations (20)-(22). Several other forms of compact support functions can be found in 
Wendland (1999) and Wu (1995). The functions in equations (26) and (27) were specifically chosen 
because they were previously adapted and used as test functions in the MLPG method (Atluri and Shen, 
2002a). 

Hybrid RBF Interpolation 

Most of the classical RBFs shown in table 1 and the compactly supported functions in equations (26) 
and (27) cannot represent simple polynomials exactly at every point in the domain of the interpolation 
(Wang and Liu, 2002a and 2002b); they can represent the polynomial values exactly only at the N scat- 
tered points where polynomial values are prescribed. Figures 2 and 3 show RBFs (the thin curves) 
obtained from the compact RBF in equation (26) with sj = 0.6 using 5 nodes in the interval — 1 < x < 1 
for a constant polynomial of value unity and a linear polynomial, respectively. The function values are 
evaluated at the five nodes and are prescribed as the uj values in equation (17). The values of u are 
evaluated by using equation (17) at 200 points in the interval — 1 < x < 1 and are plotted in these figures 
(the thick black curves). The exact functions are also shown (thick gray lines). The compact RBF inter- 
polation recovers the polynomial values exactly only at the five nodal points, and elsewhere the values of 
the polynomials deviate from the exact values somewhat significantly. The deviation from the exact 
values can, however, be reduced by increasing the number of nodes in the domain — 1 < x < 1 . 

To improve the polynomial accuracy of the solutions, Powell (1992) suggested augmenting the radial 
basis interpolation by adding polynomial basis functions as 

N m 

w(x) = Y Rj (x) aj + Y Pk( x ) z k (28) 

7=1 k=\ 



Node: 1 2 3 4 5 


Figure 2. RBF values that correspond to interpolated and exact values of a constant polynomial. 
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Node: 1 2 3 4 5 

Figure 3. RBF values that correspond to interpolated and exact values of a linear polynomial. 


-> 

x 


where Rj, aj, and N are as shown in equation (8), p(\) are polynomial basis functions, Z£ are the unkown 
coefficients associated with the kth polynomial term, and m is the order of the polynomial basis function. 
Equation (28) is written in matrix form as (Wang and Liu, 2002a) 


m(x) = R T (x) a + p T (x) z 


R T (x) p T (x) 


l a l 

I Z 


(29) 


where a and R T (x) are as in equation (12), and 

z = {zi, z 2 , Z3, .. 


J 1 


P T (x) = [p>i(x), p 2 (x), p 3 (x), ..., P m (x)\ 


(30) 


In 1-D problems, 


P T (x) = 


, 2 m — 1 

1 , X, X , X 


The interpolation is forced to pass through the N points, with the constraint 


N 


^ p k (Xj)aj = 0, * = 1,2,..., 
7=1 


m 


(31) 


imposed to guarantee unique approximation (Wang and Liu, 2002a). The set of equations to determine 
the coefficients {a} and {z} can then be written as 



where u T and Rg are defined in equations (14) and (15), and 


(32) 


P\ (X| ) 

p 2 ( x l) •' 

P ;M (X] ) 

Pl(x 2 ) 

P 2 (X 2 ) 

Pm ( x 2 ) 

Pl(x/v) 

P2( x At) •' 

Pm ( X T7) 


The unknown coefficients in equation (32) are obtained as 



The interpolating functions «(x) in equation (29) can now be rewritten as 


m(x) 


R T (x) 


P T (x) 


[G] 


-l 



N 

= X <Py( x )«/ 

7=1 

The nodal shape functions are then 


(33) 


(34) 


(35) 


N m 

9/ (x) = X R j ( X )T 77 + X Pk (x)Y (N+kYD (36) 

7=1 k = 1 

where y pq are the elements of the matrix [G] _1 . The shape functions in equation (36) possess the 
properties of equations (20) and (21). The derivatives of the shape functions are easy to evaluate using 
equation (36) as 


dtp /(x) 
dx 


dtp/(x) 

dy 


£ d^'(x) 


m 

yj / + X 


A=1 


« dK,(X| 

£1 * 


m 

yji+t 


k = 1 


dPA-(x) 

dx 


dPkW 

dy 


v (N+W) 


y{N+k)d) 


where (dRj/dx) and (dRj/dy) are as defined in equations (24) and (25). 
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( 37 ) 



Consider again the constant and linear polynomials depicted in figures 2 and 3. When a constant 
polynomial basis (p T (x) = [1]) is included in the interpolation, as described by equations (28)-(37), 
the constant polynomial of value unity (see fig. 2) is recovered exactly at all locations in the domain 
(the solution coincides with the exact solution in the figure). Similarly, when a linear polynomial basis 
(p T (x) = [1, x]) is included in the interpolation, the values of the linear polynomial (fig. 3) are recovered 
exactly at all locations in the domain. 

Interpolation Schemes — Summary 

In the previous section, two interpolation schemes were presented for C° problems: the MLS interpo- 
lation and the RBF interpolation. In this section, the FEM, MLS, and RBF interpolations are compared in 
terms of the continuity of the secondary variables. For C° problems, the secondary variable is, for 
example, ( q = du/dx). Examples of trial functions formed for each of the interpolations are presented 
schematically in figure 4. 

In the FEM, u is chosen as a piecewise linear function, as is shown in figure 4(a). (For simplicity, the 
nodal shape functions for element e are shown, and the shape functions for the remaining elements are 
omitted from the figure.) The slopes at node j, evaluated from elements e and e+ 1, are obviously 
unequal. In general, all secondary variables in the FEM are discontinuous across element boundaries 
because of the piecewise nature of the approximation for the trial function. Postprocessing techniques are 
required to achieve smooth distributions for the secondary variables. 

As discussed previously, an MLS interpolation is a scheme that fits a smooth function through an 
assumed set of fictitious nodal values. Because the trial functions formed by the MLS interpolation are 
smooth (see fig. 4(b), in which shape functions are omitted entirely from the figure), the secondary 
variables are continuous at eveiy point in the domain of the trial functions. 

In the RBF interpolation scheme, a smooth function passes through the set of actual nodal values 
obtained from the radial basis shape functions. Again, because these trial functions are smooth (see 
fig. 4(c), only shape function for node j is shown), the secondary variables are continuous. 

Beam Problem (C*) Interpolation Schemes 

This section presents the interpolation schemes used in the RPG method for beam (C 1 ) problems. The 
GMLS interpolation is discussed first; then, the shape functions for both the radial basis and hybrid 
interpolations are presented. These shape functions will be used in the beam problem implementation of 
the RPG method. 

Generalized MLS (GMLS) Interpolation 

The GMLS approximations are used for beam problems in the traditional MLPG method (Atluri, Cho, 
and Kim, 1999; Phillips and Raju, 2002; Raju and Phillips, 2002, 2003a, and 2003c). The primary 
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(a) FEM. 



(b) MLS. 



(c) RBF. 

Figure 4. Comparison of FEM, MLS, and RBF trial functions for C° problems. 
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variables for these problems are the deflection w and the slope 9 (9 = dwldx). The shape functions for w 
and 9, using a local coordinate approach (Raju and Phillips, 2993a and 2993c), are 




O) 


171 

U) = I P g (k j )[[A]" 1 [Pfft]]^ and 


g=l 


wi 

Yf 00 = I Pg&jil A] _1 [Px] T mlgj 


g = 1 


(38) 


or 




i w) = Ipg[[Ar 1 [B w ]' 


g=l 


g/ 


(39a) 


and 


m 

< ) = Ipg[[Ar 1 [B 0 ]] 

g=i 


s/ 


(39b) 


where 


[[B w ] [B e ]] = 


[P] T M [P v ] T [X] 


(49) 


In equations (38)-(49),/?(£) is a polynomial basis function, and 


[A] = [P] T MP] + [PJ T [X][P a ] 


(41) 


In equation (41), [P] and [P x ] are (n, in) matrices, and [A,] is a diagonal (n, n) matrix defined as 

-|T 


[P]= P (Si) P fe>) ••• P (£«) 


(42a) 


[Pw] = f P.v(^l) pjfe) ••• P I(5n) 


(42b) 




X = \ 


M x ) 


(43) 


K( x ) 


12 



(44a) 


where ^=Xj c -Xj ,k = 1,2, . . n. 


P T (^) = 


1 , 5 . V 


:m — 1 


and 


as 


ax 


0, 1, 2^, ... (m- \% m ~ 2 


d_ 

dx 


( 


) 


d_ 

dt, 


( ) 


(44b) 


(45) 


In C 1 problems, the MLPG method requires the first, second, and third derivatives (Atluri, Cho, and Kim, 
1999; Gu and Liu, 2001a; Phillips and Raju, 2002; and Raju and Phillips, 2003a). The first derivatives of 
\|/y are 



+ P g { [AV l [B w l x + [A]-'|BJ ) j 

§J J 


m f , , 

g=l ^ 

+ p g ( [A]-' [BqX x + [A]“ 1 [B e ] ) j 

8J J 


(46a) 


(46b) 


where [ A | J is defined in equation (7). The second and third derivatives involve considerably more 
complex expressions containing [A] -1 , [A]” 1 , [A]~\., and so on, and the detailed derivations are given 

in Phillips and Raju (2002) and are not repeated here. For thin plate problems (2-D C 1 problems), these 
derivatives become even more complicated (see Long and Atluri, 2002 for details). Evaluation of the 
derivatives is computationally expensive, which appears to be a disadvantage of using the GMLS for the 
trial functions. Thus, an alternate formulation for trial functions that is computationally simple and yet as 
accurate as the GMLS approximations is sought. The RBFs appear to be good candidates for achieving 
such a purpose because the shape functions obtained from radial basis interpolation are simpler than the 
shape functions presented previously for the GMLS. More importantly, the derivatives of the radial basis 
shape functions are simple and involve considerably fewer matrix inversion and multiplication operations 
than the derivatives of the GMLS shape functions. The RBFs are discussed next. 
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RBF Interpolation 

In the beam problems considered in this paper, the deflection w and the slope 9 = dw/dx are the pri- 
mary variables. The interpolating function for w(x) is assumed to be of the form 

w(x) = R\ (x)oj + R 2 (x)a 2 + S 2 (x)b 2 + ... +Rpj(x)a^ + S N (x)b N (47) 

where aj and bj,j = 1,2, ..., N, are unknown coefficients, Rj(x) are the RBFs, and Sj(x) = dRj(x)/dx at 
node j. Because of the direct relationship between the slope and the deflection, the approximating 
functions for 9 cannot be chosen independently from the functions for w (eq. (47)); thus, the approxima- 
tions for 9 are written as 


9(x) = 


dw(x ) 
dx 


dRi (x) dSi (x) , 
-«t + — — h 


dx 


dx 


dR 2 (x) dS 2 (x ) r 

+ a 2 + — - — b 2 + 


dx 


dx 


cIR n (x) dS N (x) 

+ a N + b N 


dx 


dx 


In matrix form, equation (47) is written as 

w(x) = Q T (x){c} 


(48) 


(49) 


where 


Q T (x) = [i?i(x) 5, (x) R 2 (x) 5 2 (x) ... R n (x) 5 w (x)] 
{c} = {«i b { a 2 b 2 ... a N b N } T 


Similarly, equation (48) is written as 


9(x) = 


dQ T (x) 

dx 


{c} 


(59) 


(51) 


where 


i7Q T (x) 
dx 

The interpolations in equations (59) and (52) are required to pass through N nodal values, which yields 

[Q B ] {c} = {d} (53) 

(2 N, 2 N) (2 TV, 1) (2 TV, 1) 


dR\ (x) dS\ (x) dR 2 (x) dS 2 (x) 


dR N (x ) dS\;(x) 


dx 


dx 


dx 


dx 


dx 


dx 


(52) 


where 


{d} T = {w 1 9 l w 2 9 2 ... w N 9^} 
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( 54 ) 



is the vector of nodal values of w and 9 at the N nodes, and 


R\ ( x \ ) 

S\ (X [ ) 

R 2 (x l) 

S 2 (x\) ... 

Rn ( x i ) 

S N (x i) 

dR\ (v| ) 

dS\ (a'| ) 

dR 2 {x i) 

dS 2 (x\ ) 

r//T v (A| ) 

dSVOq ) 

dx 

dx 

dx 

dx 

dx 

dx 

R { {x 2 ) 

S\{x 2 ) 

R 2 {x 2 ) 

S 2 (x 2 ) 

Rn ( x 2 ) 

s N ( X2 ) 

dR\ (x 2 ) 

dS\ (x 2 ) 

dR 2 (x 2 ) 

dS 2 (x 2 ) 

dR N (x 2 ) 

dS N (x 2 ) 

dx 

dx 

dx 

dx 

dx 

dx 

R\ (x N ) 

S\{x N ) 

R 2 ( x n) 

S 2 (x N ) 

Rn( x n) 

Sn( x n) 

(x, v ) 

dS\ (x N ) 

dR 2 (x N ) 

dS 2 (x N ) 

dR N (x N ) 

dS N (x N ) 

dx 

dx 

dx 

dx 

dx 

dx 


The unknown coefficients in equation (53) are obtained as 

{c} = [Q B r 1 {d} 

The interpolation for w in equation (49) can be written as 

w(x) = Q T (x)[Q B r 1 {d} 

2 N 

= X ^ ji^dj 

7=1 


where cp are the nodal shape functions 

(p(x) = Q T (x)[Q B r 1 


¥i w) W Vi 0) ( x ) V? 0) ( 


( 55 ) 


(56) 


(57) 


(58) 


From equation (58), the individual shape functions for deflection and slope, \|//' \ x ) an d V|/^(x) , are 

N 

V/ H) 0)= X (%( x )' r l(27:-l)(2/-l) + 5 , yt(x)-11(2fc)(2/-l)) 
k = 1 

(59) 

N 

V/ 9) (x)= X (^k( x ) ' ^l(2k-l)(2/) T Sk (x) • r((2fc)(2/)) 

k = 1 
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where r\,, are the elements of the matrix [Qb] 1 • The derivatives of these shape functions are easily 


Hi 

evaluated as 


_ * f dR k ( x) 

^ dx 


dx 


' r l(2/t— 1)(2/— 1) 


+ 


dS k (x) 

dx 


’ r l(2Ar)(2/-l) 


d\ ( dR k (x) 


N 


dx 


=x 

k=\ 


dx 




+ 


dS k (x) 

dx 


-\2k){2l) 


(60) 


d 2 \\i^ w \x) 
dx' 2 


= X 

k = 1 


N ( d 2 R k (x) d 2 S k (x) A 

—J r l(2fc-l)(2/-l) + — — 2 r l(2fc)(2/-l) 

t/X Ct-/V 


d~\\f^\x) 

dx 2 


N 


=x 

jfc=l 


^ d 2 R k {x) 


dx 


2 r l(2fc-l)(2/) 


+ 


d^ S k (x) 
dx 2 


T l(2fc)(2/) 


(61) 


r/ 3 t|/| w ^(x) 

dx 2 

d 3 \sp^(x) 

dx 2 


N 


= X 

k = 1 


f d 3 R k (x) d 3 S k (x ) 

~7 3 T I(2*-1X2/-1) + “ T 3 r l(2A:)(2/-l) 

CiA (tA 


tV 


ifc=l 


^ d 3 R k (x) d 3 S k (x ) 

~7 3 r l(2*-l)(2/) + — 3 r l(2*)(2/) 

ax ax 


(62) 


Hybrid RBF Interpolation 

As discussed for C° problems, to improve the polynomial accuracy of the solutions, interpolations 
involving both RBFs and polynomial basis functions are considered as 


w(x) = Q T (x) (c}+p T (x) {z} 


Q (x) p (x) 


(63) 


J z 


where {c} and Q T (x) are as in equations (50), and p T (x) are the polynomial basis functions 

p T (x) = [pi(x), p 2 (x), p 3 (x), p m (x)\ 


I i 2 m-\ 

= 1, X, X , ..., X 


and {z} are the unknown coefficients associated with p T (x) 


{z} = {zi, Z 2 , Z 3 , zj 


(64) 
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The interpolation of equation (63) is required to pass through the N points with constraints 


I 


f 

Pk(xj)c 2 j-\ 


j= l v 


+ - 


dPkixj) 

dx 


~ C 2j 


= 0 


( 65 ) 


where k =1,2,..., m, and {c} are the unknown coefficients in equations (50) (Wang and Liu, 2002a and 
2002b). Equation (65) is imposed to guarantee a unique approximation. The set of equations to 
determine the coefficients {c} and {z} is thus written as 


T TB ft 

Tb 0 J |zj 



( 66 ) 



where {d} T and [Qb] are defined in equations (54) and (55), and [Tb] is a (2 N, m ) matrix 


P] (a ) 

P2^ x \) 

••• Pm( x l) " 

dp\ (x, ) 

dp2 ( x l ) 

dp mix j) 

dx 

dx 

dx 

P\ ( x 2 ) 

P2 (. x 2 ) 

••• Pm( x 2) 

dp\ (x 2 ) 

dP2 ( x 2 ) 

dp m (x 2 ) 

dx 

dx 

dx 

Pl(x N ) 

P2( x n) 


dp\(x N ) 

dp 2 (x N ) 

dp m (x N ) 

dx 

dx 

dx 


The unknown coefficients in equation (66) are obtained as 



The interpolating function for w in equation (63) is now written as 


w(x) = 


Q T (x) 


p T (x) 



2N 

= X ( ?j( x ) d j 

7=1 


(69) 
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where (p are the nodal shape functions 


(p(x) = 


Q T « P T W 


[G] 


i-l 


V|/ ( j w) (x) y{ 0) (*) t|/ ( 2 W) (x) t|/ ( 2 0) (x) ... t| r { u\x) v|/5J } (x) 


(70) 


From equation (70), the shape functions for the deflection and slope \|/7 Vx) and t|/ 0 ^(x) may be written 
as 


N 

V/ H ) ( x ) = X ( R j ( x ) ' C(2)-l)(2/-l) + Sj ( x ) ■ C(2 i)(2/-l) ) 
7=1 


+ X Pk( x y^(2N+k)(2l-l) 
k = 1 

(71) 

N 

¥/ 0) (*) = X [ R j W ' C(2 )-l)(2/) +Sj(x)- C,(2j)(2l) ) 

7=1 


+ X Pk( x y^(2N+k)(2l) 
k = 1 

where ^ are the elements of the matrix [G] -1 . The derivatives of these shape functions are easy to 
evaluate, as in equations (60)-(62). 

MLPG Equations for Beam (C*) Problems 


Beam problems are governed by the fourth-order equation 


El -—^- = f in 0 < x < L (72) 

dx 


subjected to four boundaiy conditions, two at each end (x = 0 and x = L). The boundaiy conditions are on 
deflection w, slope 0, moment M, and shear force V, where 



dx 


M = EI 


,2 

a w 
dx 2 


V = -El 


j3 

d w 
dx 3 . 


( 73 ) 
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The essential boundaiy conditions are on w and 9, while the natural boundary conditions are on V and M. 
The boundaiy condition sets on w and V and 9 and M are disjoint, that is, if w is prescribed, V cannot be 
prescribed, and vice versa. 

The MLPG equations are derived using a weighted residual weak form of the governing equation 
(eq. (72)). The MLPG equations are (Atluri, Cho, and Kim, 1999; Phillips and Raju, 2002; and Raju and 
Phillips, 2002 and 2003a) 


j^Cnode)^ + j^(bdry)^ _ j (node) _ j(bdry) _ q 
where the superscript “bdry” denotes boundary, 

d T = {w 1 0! w 2 0 2 ... w N 0^} 


are the nodal values of deflections w and slopes 0 at the N nodes in the model of the beam, and 

j^(node) _| j^(node) 


with 
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(bdry) _ 


. (bdry) 

k ij 


k (node) _ ei 
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,(0 


dx 


dx 


J 
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dM W) d 2y(Q) 


-dx 
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-dx 


+ n x EI 


X. 


(w) 
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dx 


d\ ( f 

dx 3 

d\ ( j v) 

dx 2 
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n ( 


(0 


dx 


dx 


Xi 


(w) 


dVP 

dx 3 


X 


d 3 ur (0) 

( 0 ) X , 


J 


dx~ 


Jr 


to 


- n x EI 


d Xi 


d 2 ^^ dv {w) d 2 \|/ (0) 


dX) 


V 


dx dx 2 


dx dx 2 


ix f» dWE dx'E dWP 


dx dx 2 


dx dx 2 


<i) 


(74) 


(75a) 


(75b) 

(75c) 


(75d) 
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and 


k < M ,»_ a Kvr w 

ij _ Ixfvr xfvf L 


+ n x EI 


Y (w) Y i v (w) W 

'W A 3 ki A 3 
ax ax 


jeifV 1 .(9)‘'V/" 

,3 'w ,3 
ax ax 


<Xz{ w> ‘*< ) rfxS“’ ) ‘V, 


dx dx 

^ e) <’ 

. dx dx 


dx dx 
dx dx Jr (,) 


— n x EI 


' d 2 ^ dyff) d 2 ^f 
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q(0 
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(75 e) 


(75f) 
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where i = 1,2, ..., N and j =1,2,..., n, and n is the number of nodes in the domain of definition of the 
trial function. In these equations, %■" ) and y*f } are components of the test functions, \|//" * and V/^ 1 

are the shape functions, Q (see fig. 5(b)) is the local subdomain of the test function at node i (see 
fig. 5(a)), n x is the unit outward normal to Q ^ (n x = -1 at the left boundary point, and n x = +1 at the 
right boundary point), and are the boundary points of Q ^ (see fig. 5(b)). When coincides with 
an interior point, that point is denoted T^j, and r ( v g, f ^ . and F^y denote the boundary points 

where intersects the boundary when w, 9, M, and V are prescribed, respectively. Also in these 
equations, a w and OCg are penalty parameters to enforce the essential boundary conditions, and w, 9, 
M, and V are prescribed values of the deflection, slope, moment, and shear, respectively. See Raju and 
Phillips (2993a) for a detailed explanation of these terms. 

The system of equations presented in equations (74)-(75g) is the general set of equations valid for any 
set of trial and test functions. In this paper, a Petrov-Galerkin method is used; the test functions are 
chosen to be different from the trial functions. The choices for the trial and test functions are now briefly 
discussed. 


1 2 k 



x 


(a) ,/V-node model of beam. 



(b) Components of trial and test functions. 

Figure 5. Comparison of domains of trial and test functions. 
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Trial Functions 


The trial functions are chosen as 


N 

w(x) = ]T (vf\xywj + \\if ) (x)Q J j (76) 

./=l 

where \|/b \x) and V|/^(x) are the radial basis shape functions of equations (59), and wj and 0y are the 

nodal values of w and 0 at the N nodes (eq. (54)). Note that in MLPG algorithms that employ the moving 
least squares interpolation scheme for the trial functions, the values d (eq. (75a)) are fictitious nodal 
values d. In this paper, RBFs pass through the actual nodal values d. 

Test Functions 


The test function v is assumed as in Raju and Phillips (2003 a) as 


v(x) = l4' r) X- H ’ ) (x) + b| 0) X/ 9) W 


(77) 


where |i* H 'and (i| 0 ^ are arbitrary constants. In this paper, the %*j w \x) components of the test functions 
are chosen as power weight functions (Raju and Phillips, 2003a) 


xb’frM 



U'l 

2“ 

1- 




v s oy 



-1 4 


if 0 < dj < s a 


if dj > s a 


(78) 


with di = ||jc - Xi\\. In equation (78), s 0 is a user-defined parameter that determines the extent of the test 
functions (and hence see fig. 5). The components of the test functions chosen for 0 are the first 
derivatives of the components of the test functions chosen for the primary variable w, that is, 


X 


( 0 ) _ dli 


(w) 


dx 


(79) 


For this power function, the values of %- w \ (dy^" ] I dx). and (d~yf ] Ulx) are zero when d t = s 0 . 

As discussed in Raju and Phillips (2003a), when this test function is used, the kd node ) in equa- 
tion (75d) reduces to 


. (node) 

k ij 


1 


dh { r' d2 ^\. f 


£2l° 


dx~ 


dx ' 


-dx 


1 


£2(° 


dx 1 


dx~ 


-dx 


El 


I 


rf 2.6) ,,2 ¥ <e» 


n' 


,(0 


dx~ 


dx 1 


-dx 


J 




(0 


dx 1 


dx 1 


-dx 


(80) 


22 



Numerical Evaluations 


In this section and the sections that follow, several example problems are presented to validate the 
RPG method for beam (C 1 ) problems. First, the beam configurations and models used in the numerical 
evaluations are described, and then three patch test problems are described. Next, the compact RBF 
interpolations and their hybrid form are used in the RPG method to study the patch tests and simple 
mixed boundary value problems. The performance of noncompactly supported RBFs is presented next. 
Finally, the RPG method with a noncompact cubic function is used to study problems with complex 
loading and boundary conditions. 

Beam Configurations and Models 

A beam of constant flexural rigidity El and a length of 4/ is considered. The length 4/ is specifically 
chosen to avoid scaling by a unit length 1. Five models with 5, 9, 17, 33, and 65 nodes uniformly distrib- 
uted along the length of the beam are considered. Figure 6 shows a typical 17-node model. The distances 
between the nodes (Ax//) in these models are 1, 0.5, 0.25, 0.125, and 0.0625 for the 5-, 9-, 17-, 33-, and 
65 -node models, respectively. Numerical integration (a 12-point Gaussian integration) is used to 
integrate the terms in the system of equations shown in equations (75d) and (75f). 


|-< 4 / >~| 

I — •— 9 — 9 — 9— 9 1 — I » — t 1 — • • — • 1 

1 2 9 ->|ax|<- 16 17 




Figure 6. A 17-node model of beam. 


Patch Tests 

The radial basis MLPG (RPG) method is evaluated by applying the method to simple patch-test prob- 
lems. The problems considered are 

(a) rigid body translation: 


(b) rigid body rotation: 


w(x) = (3 0 


0 = ^ = 0 

dx 


w(x) = |3 jX 9 = Pl 


(81a) 


(81b) 


(c) constant-curvature condition: 


w(x) = P 2 x 2 /2 9 = P 2 x (81c) 

where Po, Pi, and p 2 are arbitrary constants. The third patch test is equivalent to the problem of a 
cantilever beam with a moment, M = EI{d~w/dx ~ ) = £7p 2 , applied at x = 4/. The deflection w and the 
slope 9 that correspond to problems (a), (b), and (c) are prescribed as essential boundary conditions 
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(EBCs) at x = 0 and x = 41. With these EBCs, the beam problems are analyzed by using the RPG method 
with no polynomial basis. If the RPG method recovers the exact solution at all the interior nodes and at 
every arbitrary point of the beam, then the method passes the patch test. 

Compact RBF 

The compact radial functions described by equations (26) and (27) were considered first. When using 
the compactly supported functions, equations (55)-(62) are evaluated with N = n, the number of nodes in 
the influence domain of the point x under consideration (Wang and Lu, 2002a and 2002b). This use 
of the compact functions forces the [Qb] of equation (55) to become a (2 n, 2 n) matrix that must be evalu- 
ated once for every node in the model. For these functions, the value of s a (eq. (78)) is chosen as 
(s 0 ll = 2 Ax//), and the value of sj (eqs. (26) and (27)) is chosen as (sj/l = 8A x/l). 

The RPG method with no polynomial terms was unable to reproduce the exact solutions in 

T 2 

equation (81) and thus failed the patch tests. A quadratic polynomial basis (m = 3; p (x) = [1, x, x“], 

p A (x) = [0, 1, 2x]) was used in equation (63), increasing the size of the [Qb] matrix to (2 n + m, 2 n + m). 
The RPG method with a polynomial basis (the hybrid RPG method) reproduced the solutions in equa- 
tion (81) to machine accuracy, thus passing the patch tests. 

Next, mixed boundary value problems were considered. The first problem considered was a cantilever 
beam with a tip load (fig. 7). Because the exact solution for this problem is cubic in x, the hybrid RPG 
method with a cubic polynomial basis function reproduced the exact solution. A simply supported beam 
subjected to a uniformly distributed load (fig. 8) was considered next. Because the exact solution for this 
problem is quartic in x, the hybrid RPG method with quartic basis yielded the solution exactly. 
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Figure 7. Cantilever beam with tip load. 
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Figure 8. Simply supported beam subjected to uniformly distributed load. 
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As with the finite element method and the conventional MLPG method (Phillips and Raju, 2002; and 
Raju and Phillips, 2003a), the hybrid RPG algorithm should be robust enough to yield good solutions 
when a low order polynomial basis function is used. A convergence test was conducted to study the 
performance of the hybrid method in solving the problems in figures 7 and 8. A quadratic polynomial 
basis function was used. For all models (5, 9, 17, 33, and 65 nodes), the method did not yield meaningful 
results. Thus, it was concluded that as long as the order of the polynomial basis was sufficient to repro- 
duce the solution exactly, the polynomial terms overpowered the RBFs. This condition is too restrictive, 
and hence compact radial functions are dropped from further consideration. 

Numerical Evaluations With Noncompactly Supported RBF 

The noncompact functions of table 1 are considered next. As stated previously, in these functions, 

tj=— (82) 

S J 

where rj = ||x - xj\\, and sj is some normalizing distance, usually chosen to be the entire problem domain, 
Q (here, Q is 0 < x < L ). Because sj covers the entire problem domain, [Qb] is an {IN, 2 N) matrix that 
is evaluated and inverted once. As a result, the current RPG method involves one inversion of a large 
matrix and one solution of a large system of equations ([Qb] - 1 to obtain the shape functions and a 
solution involving [K]). In contrast, the conventional MLPG method involves many inversions of smaller 
matrices to obtain the shape functions and their derivatives and the solution of a large system of equations 
that involves the [K] matrix. Because of this difference, the RPG method may not be more compu- 
tationally efficient than the MLPG method for large models. 

Upon examination of all the functions in table 1 , the cubic function (with c = 0) 

Rj(x) = t 3 j (83) 

was chosen because it is simple and gives nonzero second and third derivatives, as required by the weak 
form. (Note that Rj(x) = tj cannot be chosen because the second and higher order derivatives are zero.) 
The RPG method using equation (83) was applied to the patch tests represented by equations (81). The 
method successfully reproduced the exact solutions to machine accuracy, thus passing all the patch tests. 
Additionally, all functions of the form 


Rj(x) = tf z 1} (84) 

3 

where z > 1, performed successfully, although tj yielded the best results (Raju, Phillips, and 

Krishnamurthy, 2003b; and Raju and Phillips, 2003d). Note that when this function is chosen, the 
polynomial terms are not included. 

Next, the RPG method with the RBF in equation (83) was used to analyze mixed boundary value 
problems. In the method, a 12-point Gaussian integration was used to integrate the weak form, and the 
value of (s 0 /l), which defines the extent of the test functions (see eq. (78)), was set as {s 0 ll = 2A x/l). 
Because a noncompact function was chosen for the trial functions, the value of ( sj/l ) (eq. (82)) was set 
as (sj/l = LIT). For the cantilever beam with a tip load in figure 7, the RPG method yielded the exact 
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solution. The simply supported beam problem with a uniformly distributed load (fig. 8) was analyzed by 
using 17-, 33-, and 65-node models. The maximum deflection values (i.e., the deflection at (x = i/2)), for 
these three models, obtained by using the RPG method and the conventional MLPG method with a quad- 
ratic polynomial basis function, are compared in table 2. In the MLPG method, a 20-point Gaussian 
integration was used to integrate the weak form (Raju and Phillips, 2003a), the value of (s 0 /l) was set as 
(s 0 /l = 2A x/l), and the value of (sj/l) was set as (sj/l = 8A x/l). From table 2, it is seen that the RPG 
method yields solutions that are as accurate as the solutions obtained from the conventional MLPG 
method. 


Table 2. Comparison of Maximum Deflection 
Values Obtained Using RPG and MLPG 
Methods (ii'Exact = 5qL 4 /3S4EI) 


Model 

(w/wExact) at x = L/2 

RPG 

MLPG 

17-node 

0.982 

0.993 

33-node 

1.002 

1.012 

65-node 

1.003 

1.015 


The results obtained for deflection, slope, moment, and shear by using the 65 -node model are pre- 
sented in figure 9. In this figure, the RPG results are compared to the exact solution and to the solution 
obtained when using the conventional MLPG method with a quadratic polynomial basis function. The 
RPG values for deflection, slope, and moment are as accurate as the MLPG values and are in excellent 
agreement with the exact values. In addition, the RPG values for shear converge with model refinement. 
The MLPG solution for the shear is erratic and is not shown in figure 9. In the MLPG solution, the 
quadratic basis function is insufficient to accurately calculate the third derivatives for this problem, and 
the method cannot recover the values with model refinement; the solution for the shear converges only as 
the order of the basis function is increased to quartic (Phillips and Raju, 2002). Similar results were 
obtained by using the 17- and 33-node models; these results are not shown. The results discussed for this 
problem verify a perceived advantage of the RPG method over the MLPG method; namely, accurate 
solutions are obtained with simpler evaluation of the shape function derivatives. 

Problems With Complex Loading and Boundary Conditions 

The RPG method with the RBF in equation (83) was next applied to the complex beam problems with 
the loading and boundaiy conditions shown in figure 10. The problems considered were (a) a cantilever 
beam (L = 21) long with a uniformly distributed load on a portion of the beam shown in figure 10(a); (b) a 
continuous beam (L = 41) long with the additional support in the center of the beam (fig. 10(b)); and (c) a 
fixed continuous beam ( L = 31) long with a hinge, as shown in figure 10(c). 

The RPG solution (with (s 0 /l = 4 Ax//)) for the cantilever beam problem (problem (a)) exhibited con- 
vergence with model refinement. These results are consistent with those reported in Raju and Phillips 
(2003a), in which this problem was studied by using the conventional MLPG method. The exact MLPG 
and RPG values for deflection and moment for this problem were obtained by using a 65-node model and 
are compared in figure 1 1 (with parameters as reported previously for the simply supported beam prob- 
lem). The RPG method handled the load discontinuity well and yielded results in overall agreement with 
the exact solutions. 
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(b) Secondary variables (Note: Fmlpg not shown). 


Figure 9. RPG, MLPG, and exact solutions obtained using 65-node model for simply supported beam 
subjected to uniformly distributed load. 


For the continuous beam problem (fig. 10(b)), the RPG method yielded very accurate results for both 
the primary and secondary variables and handled the discontinuity caused by the additional support well. 
The MLPG, RPG, and exact moment distributions are shown in figure 12. For this problem, the MLPG 
method required a large number of nodes to obtain an accurate solution (Raju and Phillips, 2003a), but the 
RPG method yielded accurate results when the smaller nodal (17-node and 33-node) models were used. 
These results again verify perceived advantages of the RPG method over the MLPG method; accurate 
results can be obtained by the RPG method for comparable computing efforts. 

Next, the hinge problem shown in figure 10(c) was considered. Because the hinge at (x = /) (see 
fig. 10(c)) cannot admit the moment, there is a slope discontinuity at (x = I). Thus, there are two separate 
slopes, 0 L and 9 R , just to the left (x = l~) and just to the right (x = l + ) of the hinge, respectively, while the 
deflections, w L and w R , are equal. Three models — 14-node, 26-node, and 50-node — were considered. 
The 14-node model consisted of 5 nodes in the left segment (0 <x</) of the beam and 9 nodes in the 
right segment (/<x<3/) of the beam. Similarly, the 26-node and 50-node models were composed 
of 9 nodes in the left and 17 nodes in the right, and 17 nodes in the left and 33 nodes in the right, 
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(a) Cantilever beam with uniformly distributed load on portion of beam. 
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(b) Continuous beam subjected to uniformly distributed load. 
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(c) Fixed continuous beam with hinge. 


Figure 10. Cantilever and continuous beams with complex loading and boundary conditions. 

respectively. The problem was analyzed by using the MLPG and RPG methods. For both methods, 
the value of (s 0 /l) was set as ( s a /l = 2A x/l). For the MLPG method, the value of (sj/l) was set as 
(sj/l = 8 Ax//), and for the RPG method, the value of (sj/l) was set as (sj/l =1) for the left segment and 
(sj/l = 2) for the right segment. 

The RPG and the exact solutions for the deflection and slope at x = /, 21, and 31 are compared in 
table 3. In this table, the RPG solutions that were obtained with the three nodal refinements are pre- 
sented. Accurate solutions are obtained with the coarse model, and the solutions improve with nodal 
refinement. 
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Figure 11. RPG, MLPG, and exact solutions obtained using 65-node model for cantilever beam with 
uniformly distributed load on portion of beam. 



Figure 12. Comparison of RPG (33-node model), MLPG (65-node model), and exact solutions for moment of 
continuous beam subjected to uniformly distributed load. 


Table 3. Comparison of RPG Results and Exact Solutions at Several Locations Along Beam 


Model 

w/w Exact (x = l) 

0 L/0 Exact(' t = O 

0 i? /0ExactU = O 

0 /0 Exact(^ = 2/ ) 

Q/^ExactO^ 3/ ) 

14-node 

0.9769 

0.9748 

0.9679 

0.9774 

0.8944 

26-node 

0.9814 

0.9789 

0.9798 

0.9817 

0.9656 

50-node 

0.9880 

0.9867 

0.9868 

0.9882 

0.9755 
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Figure 13. RPG, MLPG, and exact solutions for deflection of fixed continuous, hinged beam subjected 
to uniformly distributed load. 



Figure 14. RPG, MLPG, and exact solutions for moment distribution of fixed continuous, hinged beam 
subjected to uniformly distributed load. 


The RPG and MLPG solutions obtained with the 50-node model are compared in figures 13, 14, and 
15 for deflection, moment, and shear force, respectively. For the deflection and moment distributions, 
both the RPG and MLPG solutions agree very well with the exact solutions. The MLPG method is 
unable to accurately represent the shear force; hence, this solution is not included in the figure 15 com- 
parison. The RPG solution, on the other hand, agreed very well with the exact solution. 
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Figure 15. RPG and exact solutions for shear force distribution of fixed continuous, hinged beam subjected 
to uniformly distributed load (Note: MLPG not shown). 

Concluding Remarks 

A radial basis function (RBF) implementation of the meshless local Petrov-Galerkin (MLPG) method 
in the study of Euler-Bernoulli beam problems was presented. Like the conventional MLPG method, this 
radial basis variation (RPG) is based on the local weak form developed from the classical weighted 
residual form of the fourth-order governing differential equation. In this method, rather than generalized 
moving least squares (GMLS) interpolations, RBFs were used to develop the trial functions, and test 
functions were chosen as simple weight functions, as in the conventional MLPG method. Two types of 
RBFs, compact and noncompact, were considered. Compact RBF interpolations were developed with and 
without including polynomial basis function terms, while noncompact RBF interpolations without 
polynomial terms were considered. 

The compactly supported RBFs did not perform well without polynomial terms in the computations. 
When polynomial terms were included, the compactly supported RPG method passed the patch tests. 
However, the method did not yield meaningful results for mixed boundary value problems unless the 
order of the polynomial basis function was of the same order as the exact solution of the problem. This 
result restricts the use of the method. Also, polynomial terms tended to overpower the RBFs; hence, the 
use of compactly supported RBFs presented in this paper is not recommended in the RPG method. 

The noncompactly supported cubic RBF performed very well when no polynomial terms were in- 
cluded. The use of RBFs rather than the traditional moving least squares interpolations reduced the 
computing effort required to solve problems; substantially fewer matrix inversion and multiplication 
operations were required by the RPG to evaluate the derivatives of the shape functions. When noncom- 
pactly supported RBFs were used, the RPG method involved one inversion of a large matrix to evaluate 
the shape functions. The RPG method with the cubic RBF yielded accurate solutions for the problems 
studied; the method passed all the patch tests, and for mixed boundary value problems, achieved the same 
accuracy in results as the conventional MLPG method, with comparable computing effort. The RPG 
method with a cubic radial basis also yielded very good results in contrast to the conventional MLPG 
method, which could not yield consistently good results for problems with complex loading and boundary 
conditions. The accuracy of solutions obtained by the RPG method, combined with the computational 
simplicity of using the RBFs rather than the GMLS interpolations to approximate the trial functions, 
makes the RPG method a very attractive variation of the MLPG method. 
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